home *** CD-ROM | disk | FTP | other *** search
/ Mac Magazin/MacEasy 21 / Mac Magazin and MacEasy Magazine CD - Issue 21.iso / Wissenschaft & Technik / yorick12vr1-nofpu folder / include / gcd.i < prev    next >
Text File  |  1995-12-12  |  4KB  |  142 lines
  1. /*
  2.    GCD.I
  3.    GCD, LCM, and prime factorization routines.
  4.  
  5.    $Id$
  6.  */
  7. /*    Copyright (c) 1995.  The Regents of the University of California.
  8.                     All rights reserved.  */
  9.  
  10. func gcd(a, b)
  11. /* DOCUMENT gcd(a,b)
  12.      returns the GCD (greatest common divisor) of A and B, which must
  13.      be one of the integer data types.  A and B may be conformable
  14.      arrays; the semantics of the gcd call are the same as any other
  15.      binary operation.  Uses Euclid's celebrated algorithm.
  16.      The absolute values of A and B are taken before the operation
  17.      commences; if either A or B is 0, the return value will be 0.
  18.    SEE ALSO: lcm, is_prime, factorize
  19.  */
  20. {
  21.   a= abs(a);
  22.   b= abs(b);
  23.   c= min(a, b);
  24.   a= max(a, b);
  25.   b= c;          /* simplifies c=0 case */
  26.  
  27.   if (dimsof(a)(1)) {
  28.     /* array case */
  29.     for (list=where(c) ; numberof(list) ; list=where(c)) {
  30.       b(list)= bl= c(list);
  31.       c(list)= a(list) % bl;
  32.       a(list)= bl;
  33.     }
  34.  
  35.   } else {
  36.     /* scalar case can be less baroque */
  37.     while (c) {
  38.       b= c;
  39.       c= a % b;
  40.       a= b;
  41.     }
  42.   }
  43.  
  44.   return b;
  45. }
  46.  
  47. func lcm(a, b)
  48. /* DOCUMENT lcm(a,b)
  49.      returns the LCM (least common multiple) of A and B, which must
  50.      be one of the integer data types.  A and B may be conformable
  51.      arrays; the semantics of the lcm call are the same as any other
  52.      binary operation.
  53.      The absolute values of A and B are taken before the operation
  54.      commences; if either A or B is 0, the return value will be 0.
  55.    SEE ALSO: gcd, is_prime, factorize
  56.  */
  57. {
  58.   d= gcd(a, b);
  59.   /* two potential problems: zero divide and overflow - handle the
  60.      first but not the second */
  61.   return abs(a*b)/(d+!d);
  62. }
  63.  
  64. func is_prime(x)
  65. /* DOCUMENT is_prime(x)
  66.      return non-zero if and only if X (which must be a scalar integer)
  67.      is prime.  May return a false positive if X is greater than about
  68.      3e9, since at most 20000 candidate factors are checked.
  69.      The absolute value of X is taken first; zero is not prime, but 1 is.
  70.    SEE ALSO: gcd, lcm, factorize
  71.  */
  72. {
  73.   x= long(abs(x));
  74.   if (x<2) return x==1;
  75.   /* make a list of factors which includes 2, 3, and all larger
  76.      odd numbers not divisible by 3 less or equal to sqrt(x) */
  77.   top= min(long((sqrt(x)+8.)/3.+0.5), 20000);
  78.   factors= ((3*indgen(2:top))/2)*2 - 7;
  79.   factors(1:2)= [2,3];
  80.   return allof(x%factors);
  81. }
  82.  
  83. func factorize(x)
  84. /* DOCUMENT factorize(x)
  85.      return list of prime factors of X and their powers as an n-by-2
  86.      array.  May include a large non-prime factor if X exceeds 3e9.
  87.      In any event, product(result(,1)^result(,2)) will equal abs(X).
  88.      X must be a scalar integer type.
  89.    SEE ALSO: gcd, lcm, is_prime
  90.  */
  91. {
  92.   x= long(abs(x));
  93.   if (x<2) return [[x],[1]];
  94.  
  95.   /* first get rid of any prime factors less than 102 */
  96.   primes= [2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,
  97.        71,73,79,83,89,97,101];
  98.   primes= primes(where(!(x%primes)));
  99.   powers= _fact_extract(primes, x);   /* returns "deflated" x */
  100.  
  101.   if (x>1) {
  102.     /* large prime factors require a less direct approach */
  103.     top= min(long((sqrt(x)+2.)/3.+0.5), 20000);
  104.     if (top>=35) {
  105.       /* trial divisors are all odd numbers 103 or greater which are
  106.      not divisible by three and do not exceed sqrt(x) */
  107.       trial= ((3*indgen(35:top))/2)*2 - 1;
  108.       /* discard all trial divisors which do not divide x */
  109.       trial= trial(where(!(x%trial)));
  110.       /* the smallest remaining divisor must be prime - remove it and
  111.      all its subsequent multiples to find the next prime divisor
  112.      and so on until the list contains only primes */
  113.       list= trial;
  114.       for (n=0 ; numberof(trial) ; ++n) {
  115.     list(n+1)= trial(1);
  116.     trial= trial(where(trial%trial(1)));
  117.       }
  118.       if (n) {
  119.     trial= list(1:n);
  120.     grow, primes, trial;
  121.     grow, powers, _fact_extract(trial,x);
  122.       }
  123.     }
  124.     if (x>1) {
  125.       grow, primes, [x];
  126.       grow, powers, [1];
  127.     }
  128.   }
  129.  
  130.   return [primes, powers];
  131. }
  132.  
  133. func _fact_extract(primes, &x)
  134. {
  135.   if (is_void(primes)) return [];
  136.   /* first get largest power of each prime less than or equal x */
  137.   powers= long(log(x)/log(primes)+1.e-6);
  138.   factors= gcd(primes^powers, x);
  139.   x/= long(exp(sum(log(factors)))+0.5);
  140.   return long(log(factors)/log(primes)+0.5);
  141. }
  142.